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Abstract. We examine the modes admitted by the Mestel disk, 
a disk with a globally flat rotation curve. In contrast to previous 
analyses of this problem by Zang ( 1976| ) and Evans & Read 
(1998a, 1998bi, we approximate the orbits to obtain almost 



closed expressions for the kernel of the integral equation gov- 
erning the behaviour of the modes. Otherwise we, like them, 
follow Kalnajs' programme to simultaneously solve the Boltz- 
mann and Poisson equations. 

We investigate the modes admitted by both the self-con- 
sistent and a cut-out Mestel disk, the difference being that in 
the latter, a part of the matter in the disk is immobilised. This 
breaks the self-similarity and produces a pronouncedly differ- 
ent picture, both technically and in terms of the disk properties. 
The self-consistent disk is governed by a Cauchy integral equa- 
tion, the cut-out disk by an integral equation that can be treated 
as a Fredholm equation of the second kind. 

In general, our approximation reproduces the results of the 
previous works remarkably well, yielding quantities mostly 
within 5% of the values reported by Zang and Evans & Read 
and thus also the basic result that in a "standard" cut-out disk, 
only one-armed modes are unstable at the limit of axisym- 
metric stability. In the self-consistent disk, relatively compact 
expressions for the kernel allow an intuitive understanding of 
most of the properties of neutral (non-rotating, non-growing) 
modes there. We finally show that self-consistent Mestel disks 
do not admit growing or rotating modes in this sort of stellar- 
dynamical analysis. 

Key words: Galaxies: spiral - Galaxies: structure - stellar dy- 
namics 



1. Introduction 



It was not long after Kalnajs (1971, 1977, and references 



therein) and others had laid out a programme to perform stabil- 
ity analyses for galactic disks that it was applied to the Mestel 



disk by Zang (1976; henceforth cited as Zang) in a Ph.D. The- 



sis supervised by Alar Toomre. This disk, characterised by 
an infinitesimal height and a globally constant rotation curve, 
lends itself as a model for disks that show spiral structure 
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mainly for two reasons. For one, its structure is so simple that 
one of its distribution functions has made it into a textbook 



example (Binney & Tremaine 1987); secondly, and probably 



more importantly, flat rotation curves seem to be rather preva- 
lent among ( high su rface brightness) spiral galaxies (e.g., Per- 
sic & Salucci 1991 ), and it has long been suspected that flat ro- 
tation curves and strong, "modal-looking" spiral structure are 



somehow linked (e.g., Biviano et al. 1991) 



One of the drawbacks of the Mestel disk is that the orbits in 
it cannot be expressed in elementary functions. Zang circum- 
vented this problem essentially through numerical integration 
and some elegant devices exploiting the self-similarity of the 
Mestel disk (i.e., orbits of a given eccentricity look like scaled 
copies of each other regardless of the angular momentum of 
a particle on them) to curb down on the computational effort. 
While the self-similarity is welcome in simplifying the calcu- 
lation, it also implies that a self-consistent Mestel disk does not 
possess any scales. In particular, the absence of a time scale al- 
ready leads to the expectation that discrete modes will be hard 
to construct. 

Zang found that the Mestel disk is a rather bewildering con- 
struct. His formalism leads to a singular integral equation that 
is quite difficult to treat numerically. However, the structure 
of its kernel as well as physical arguments suggest that there 
are no discrete modes apart from rather exotic non-rotating and 
non-growing ones, implying that either no growing modes ex- 
ist or modes become unstable regardless of their growth rate 
and pattern speed. Thus, his investigations concentrated on a 
variant of the Mestel disk in which a certain part of the disk 
is declared "immobile", the cut-out disk. This tampered sys- 
tem is furnished with a length scale, admits discrete modes and 
becomes treatable with his numeric technique. 

But even the cut-out disk behaves contrary to the precon- 
ceptions of the mid-70ies. While Zang would have liked to see 
a mildly unstable two-armed mode, typically the only mode 
that is growing in an axisymmetrically stable disk is a one- 
armed mode, and for this, a stability limit cannot be found. 
Even when he made the inner cut-out so steep that m = 2 
modes emerg ed in axisymmetrically stable disks in a way 
that Toomre ( 1977 ) describes as unrealistic, these one-armed 
modes, again in Toomre's words continued to "plague" him. 
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Despite these somewhat unexpected findings and the fact 
that is was only fragmentarily published in Toomre (1977) and 



Toomre ( 198 1 1, Zang's thesis became a frequently-cited classic 



and the Mestel disk developed a life of its own. Recently, Evans 



& Read (1998a, 1998a; see also Read (1997; henceforth cited 



as ERI and ERII) generalised Zang's work to the entire family 
of physically plausible disks with surface densities following 
a power law, and Goodman & Evans ( 1999| ) once more tack- 
led the issue of rotating modes the untampered (self-consistent) 
Mestel disk. 

This work set out as an attempt to avoid the numerical inte- 
grations that had to be employed by Zang and ERI by approxi- 
mating the orbits and eventually deriving a closed kernel of the 
integral equation governing the problem. Among the merits of 
the availability of such a kernel is that an analytic expression 
for the kernel might allow a decision whether or not the self- 
consistent Mestel disk admits rotating modes. Also, it might 
facilitate an analysis of the processes at work at various res- 
onances of the modes in the cut-out disk in the context of a 
full disk without having to take recourse to somewhat daring 
techniques to extend WKBJ approximations across resonances 
(e.g., Mark [1976|). 



The organisation of this paper is as follows: In section 2, 
we derive the orbits and the Hamiltonian of the disk. Section 3 
basically follows Kalnajs' programme that will in the end yield 
an (almost) closed expression for the kernel of an integral equa- 
tion governing the modes in the disk. To assess the validity of 
this kernel, we will compare the results reported in ERII with 
ours in section 4. Finally, we revisit the issue of modes in the 
self-consistent disk in section 5. 



2. The Orbits 

The equations of motion in the Mestel disk are 



f = rO" — 1/r 
L, = = const. 



(1) 



Here, r and 9 denote polar coordinates, and Lz is the z- 
component of the angular momentum per unit mass. Units have 
been chosen such that the circular velocity Vc = 27rGSo = 1, 
where G is the gravitational constant, and Eg is the surface 
density at a reference radius that can be set to unity because of 
the self-similarity of the disk. 

We seek closed expressions for the solution of the equations 
([!]). Since the exact solutions cannot be expressed in elementary 
functions, we need to employ some approximation. A straight- 
forward linearisation in configuration space - the epicyclic ap- 
proximation - fails for the Mestel disk, since the Hamilto- 
nian derived from it does not contain contributions from non- 
circular motion components. Therefore we linearise in Fourier 
space, a procedure generally known as the method of harmonic 
balance (Bogoliubov & Mitropolski 1961 ). 

Inserting one of the equations (111) into the other yields 



Ll/r^ - 1/r. 



We now substitute r{t) = i?o(l + Ccos(k<)), making i?o the 
reference radius of a given orbit, C a measure for its deviation 
from circularity, and k an epicyclic frequency. After this, the 
first Fourier coefficients of (0) with respect to t are 



Li{i+e/2)-Ri{i~ef 



5/2 



ai = - 
fei = 



2Ri{i-ey{V^-^)+^Lie 



V 5/2 



(3) 

(4) 
(5) 



where and hi denote the cosine and sine coefficients, respec- 
tively. 

Dropping all higher Fourier coefficients and substituting 
the result back into (g), one finds by comparing the terms inde- 
pendent of t 



and by comparing terms proportional to cos{Kt) 



V2 

CRo' 



1-2 



e + 2 



(6) 



(7) 



The circular frequency of the guiding centre fl results from a 
Fourier transformation of the second equation in (|l]) with (^ 
inserted and is given by 



3/2 ■ 



(8) 



We can now derive the Hamiltonian for a mass point in 
the Mestel disk in action-angle coordinates by simply insert- 
ing (^), (^, and (^) into an Hamiltonian written in Cartesian 
coordinates. Assuming orbits not too far from circularity, one 
can linearise with respect to (. After again dropping all Fourier 
components of order 2 or higher, one is left with 



(9) 



where J = k('^Rq/2 is the energy in the epicyclic motion di- 
vided by K and thus the second integral. This clearly is the sim- 
plest Hamiltonian that could have anything to do with a Mestel 
disk and indeed the Hamiltonian that would result from a clas- 
sic epicyclic approximation if it worked. As a matter of fact, at 
this point one has k = V2/ and ft = 1 /Rq, reproducing the 
result of an epicyclic analysis. This Hamiltonian has already 
been used by Collett et. al. ( |1997| ). 

The set of approximated action-angle coordinates Lz, J, 
wi :~ Kt,W2 ■— is correct to first order in (, and the canoni- 
cal equations are correct to first order in ( as well. As is demon- 



(2) 



strated in Demleitner (2000), the Jacobian of the transforma- 
tion between Cartesian coordinates and this set is 1 + OiC"^), 
showing that the transformation is indeed canonical to the re- 
quired order 
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Because of the two approximations (only keeping the low- 
est Fourier coefficients and linearising with respect to Q, it is 
not a priori clear that the disk described by this transforma- 
tion and the resulting Hamiltonian (Q) does indeed behave like 
a Mestel disk. It is one of the purposes of this work to show a 
posteriori, at is were, that this is the case. 

3. Derivation of the matrix equation 

With the Hamiltonian and the transformation equations, we can 
now follow the programme outlined in the series of papers by 
Kalnajs cited above. The first step is to establish a base for the 
potential perturbation. We choose 



a cos(toi) — \/2m sin(t«i))y 



(10) 



Here, m is the number of circumferential maxima of the spiral 
pattern (the number of arms), and a parametrises how closely 
the spiral is wound, where a spiral with larger a is more trailing 
provided that ui > 0. This expression basically is the well- 
known logarithmic spiral e'" '"('')+'™^^"^* transformed into 
our action-angle coordinates with harmonics in the angle co- 

— 1/2 

ordinates higher than one disregarded. The factor Lz was 
added to ensure definiteness in the integration of the angular 
momentum. This analogy motivates following ERI in naming 
a the logarithmic wave number and rn the angular harmonic 
number We assume that all potential perturbations can be rep- 
resented as expansions of the form 



E 



(uii, W2, J, Lz) da. 



(11) 



The next step in Kalnajs' programme is to solve the lin 
earised Boltzmann equation 

dfi 



dt 



(12) 



written with Poisson brackets [ . , . ] for an elementary pertur- 
bations (pa^miwi, W2, J, Lz) of the form (p^. Here, /i is the 
unknown perturbation of the distribution function, and /q an 
equilibrium distribution function; for /o, we use the distribu- 
tion function given by Binney & Tremaine (1987i, which, in 
our units and coordinates, and after applying Stirling's formula, 
is 



k{Lz,J) 



F 



F 



■ exp 



\/2J~ 
oHTz 



with 



Expanding the Poisson brackets in ( [l2| ) yields 



(13) 



dt Lz dwi Lz dw2 



o^lT 



{X(w\) — ma^^ exp iX{w\) 



V2J 



(14) 



with an auxiliary function X{wi) — \J \/2J / Lz {a cos(ti;i) — 

\/2m sin(wi) j . On the right-hand side of this equation, a term 

of order has been neglected. 

To solve this equation, we perform a separation of variables 

by 

/iK,«;2,t)=.9iK)e™'"^-''^* (15) 

and arrive at an inhomogeneous ordinary differential equation, 
the homogeneous solution of which is 



i{LjLz—m)'W\ I y/2 



(16) 



The general solution can be obtained via a variation of con- 
stants, where the integration constant is fixed by the condition 
that the distribution function must be 27r-periodic, gi{wi) = 
gi{wi + 2tt). 

Before writing down the result, we note that in ( p^ the 
terms with in the exponent and its coefficient are quite sim- 
ilar to each other. This invites a partial integration in the result- 
ing expression, which after some algebra leads to 




i exp[iX{wi)) 



/ exp(iX{wi + w'l) — irjw'i) dw[ 



(17) 



In this expression, we have abbreviated ?] = {ojLz — m) /\/2. 

Substituting ( [T7| ) into ( [T^ yields the perturbation of the dis- 
tribution function under an elementary potential perturbation 
4'a,m- Since (|l2|) is linear in 4>a,m,, the response to an arbi- 
trary perturbation <J)™p°^'='1 of the form dill ) can immediately 
be computed by summing over all a and m. Thus, the solution 
of (na) can be written as 



response 



imposed 



(18) 



with a linear integral operator C. 

Further pursuing Kalnajs' programme, we attempt to find 
a linear combination of elementary perturbations such that 
jicsponso ^^^ggg jjjg (jjimposcd started with. To do this, the 

Poisson equation has to be solved simultaneously to the Boltz- 
mann equation. Since £ and the Laplacian occurring in the 
Poisson equation do not commute, their eigenfunctions are dif- 
ferent, and a straightforward attempt to solve the two equations 
simultaneously will result in an integro-differential equation 
that will be very hard to solve indeed. 

To avoid this, one converts the system of equations to a 
matrix equation. The idea here is to compute potential-density 
pairs, i.e., two sets of mutually orthogonal basis functions (pa.m 
and ^a,rn for each of the perturbations in density and poten- 
tial where each pair (j>a,m, Ma m satisfies the Poisson equation. 
When one expands the perturbations in the potential and the 
surface density in these functions, the Poisson equation simply 
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reads <!>,„ (a) = Mm (a) in the coefficients $,„ (a) and Aim (a) 
of these expansions. 

The basis (/>q,,„i is already fixed by (p^, and the correspond- 
ing density base function can be found using the scheme of 



Clutton-Brock (1972). For our calculation, however, the only 



important property of the potential-density pair the scalar prod- 



uct of c 



and pia,m- This scalar product is of course invariant 



under canonical transformations, and to the the required order 
Clutton-Brock's result 



de / r dr (jf^ mP-a' 
Jo 



ra.m' 



for logarithmic spirals holds. Here, 



Km{a) 



|r(! 


+ 




1 • \|2 


ttG r( 


1 


-f \m 


1 ■ \ |2 

+ 2«")l 



(20) 



is a quantity closely related to ERI's Kalnajs factor. 

After expanding the potential and surface density perturba- 
tions in (p^, one can project both sides on (j)a' ,m'- Exploiting 
(19) and inserting the (now trivial) Poisson equation, (|l8|) takes 
the form 



A 2 / /\ ^response/ /\ 

47r Km' {a )$,„,^ [a ) 



da I dLz / dJ 

-oo Jo Jo 



dwi 



oo 

E 

m— 



^^^^imposcd(^)(^^^^^^^*^ (21) 



Equating ^response ^j^^ (jjimposed^ jutggral equation results 
the solutions of which, if they exist, are self-consistent modes. 

The integrand of ( pH ) can be written as a sum Ii + /2 of a 
part without a further integral. 



/i = $ 



J-iz i(m-m')w2 

e '^-^^ 



Air^a'^ Km' ja') 

(^(a — a') Q,os{w\)---\pl{ri\ — '}Ti^ sin(i(;i)^ 



(22) 



and one in which a further integration over the auxiliary angle 
w\ has to be performed, 

/■'"./■^ ^/2FLL^"-"')-'(ma2 + V2?7) 

h = / dw\l^a^r^ 



-0 87r2CT2i^„'(a')(e-2'^^'J - 1) 

i^J cos(ii;^ +^1 ) — a' cos(wi) — y^ m sin(w[ -}-wi)^ 



(23) 



The quadratures over W2, wi, and J, carried out in that or- 
der, are rather straightforward. Details can be found in Dem- 
leitner (2000). When integrating over w[ one has to evaluate 
an expression of the type 



2tt 



(24) 



which can be done by expanding the second factor of the inte- 
grand into modified Bessel functions Ik using 



-.cos(.-K) ^ + 2 V cos{k{v - w[)) 



k=l 



(25) 



3rr- 
2ti- 



: ^> 

In (1 - ^/2i) In (1 + ^() 



Fig. 1. Poles of the integrand of (^6|) and a sketch of the integra- 
tion path. The diamonds in the left column mark the positions 
of the poles at ^{u) = ln(l - \/2^) in the case ^ > 1/V2. 



(Abramowitz & Stegun 1972). Unfortunately, the series intro- 
duced by this operation does not in general collapse and will 
be present in the kernel. 

The integration over the angular momentum is the most in- 
volved one, resembling the calculations to compute the angular 
momentum function in Zang's formalism. For the summands 
of the series resulting from the w'l integration, one has to treat 
integrals of the structure 



g(i:7+A+l)« 



2 (2^2 - (1 - e«)2; 

where 



■ du, 



(26) 



(27) 



(note that this substitution only works for nonzero m and to), 
7 = a — a' and ^ = k /m. The integer fc is a summation index 
varying between and oo, and the complete integral requires 
this expression to be evaluated for A e { — 1, 0, 1}. 

We will compute this integral using the residue theorem, 
and thus the first step is to find the zeroes of the denominator 
of the integrand. These are at 



,(1.2) 



= ln(l±\/20 + 27rij 



(28) 



for integer j. These zeroes are well-defined, since ^, as the ra- 
tio of two integers, can never be equal to l/\/2- Let for the 
moment ^ < 1/ -\/2 and 7 > 0. 

There are three more or less technical problems in the ap- 
plication of the residue theorem to this integral. The first is that 
at least one pole is on the real axis and hence on the integration 
contour The tempting idea to take the (existing) Cauchy princi- 



pal value leads to erroneous results. As Landau (1946) noted in 
a plasma-physical context, the correct reasoning is that if a sys- 
tem supports a non-growing perturbation at present and started 
in equilibrium, the perturbation must have been a growing one 
at some time in the past. Since it has evolved continuously, the 
integration contours of the problem must have been continu- 
ously deformed as well and thus cannot have jumped across 



M. Demleitner & B. Fuchs: A new approach to the problem of modes in the Mestel disk 



5 



poles. A growing perturbation corresponds to an lo with a posi- 
tive imaginary part, which with the substitution ( p7[ ) means that 
the integration contour moves into the upper half-plane. There- 
fore, as 0, one has to leave the zeroes on the real axis 
out of the contour, as is sketched in Fig. ^ 

The second problem is that as the size of the integration 
path R ^ oo (cf. Fig. |l]) one cannot continuously deform the 
contour since the poles are in the way, and the same poles make 
it hard to establish that the contribution of the line at 3u — iR 
vanishes. The standard technique to avoid this is to write down 
the series of residues and show that integral over the difference 
between this series and the integrand vanishes as i? ^ cx). 
Since this differences has no poles, the contour can be readily 
deformed. 

The last problem is that for A = ±1 the contributions at 
diu — ±R do not necessarily vanish. Indeed, one finds that 
these "border terms" 




X = -l,u 
A = 1, u ^ 



(29) 



Of course, the limits occurring in these expressions do not 
exist. However, it should be remembered that this integral is 
being evaluated as part of the kernel of an integral equation. 
Therefore, the limits can (and have to) be taken in the sense of 
a distribution limit, lima;^±oo e^'^^ /j — ±i7r(5(7). 

With these tools, one can go about evaluating the residues 
of the integrand, which can easily be done using res(/, z) = 
lim„_^2 f{u){u — z) and an application of de FHospital's rule. 
We will give the end result below. 

A few words about the other cases are in order. First, 
if ^ > l/\/2, the poles at u^p move away from the real 
axis (cf. Fig. [l|). One can keep the expression resulting for 
^ < l/\/2 by choosing an appropriate branch of the log- 
arithm in the result. Basically, one has to see to it that the 
u'p in Fig. |l] "shift down" rather than "up" by it:. To en- 
sure this, one sets ln(— 1) = —iir, and with ln(l) = a nat- 
ural branch cut is along the positive imaginary axis, so that 
f > ln(e*"^) > -37ri/2. 

If 7 < 0, one has to integrate through the lower half-plane, 
and consequently the poles on the real axis are added into the 
residue sum. The resulting expression is, of course, identical to 
the one for 7 > 0. 

After carrying out all integrations and the summation of m 
(which is trivial because the integration over W2 contributes a 
(5„i_,„', thus decoupling different angular harmonics), ( plj ) takes 
the form 

1*00 

=-^(a')*m'(«')+ / /C(a,a')$™'(a)da. (30) 



The summand with T originates from the terms with 5- 
functions and computes to 



E3l^l4?^("'^ + 2.n'^)))). (31) 



— 2fc2 
fc=i 

The kernel itself becomes 



JC (a, a) = —i 



1-1/0 



20-2gCr2(a2+Q'2+4m'2)/4+i ^^27r(Q-Q') _ 



r(i + im+i*a')| 



m ^ '(J lo(t) 



|r (I + im+ 

OC y 

- ^ Ife(r) ( {V2k - mV2)A^e('("-"')-i) H"^'-V2k) 



k=l 



-iV2k + mV)A-'=e('("-"')-i) M™'+v^fe) 



(32) 



where we have abbreviated 



A = 



y + 2TO'2)(a'2 + 2m'2) 

aa' + \/2TO'(\/2m' + i{a - a')) 
^(2m'2 + a2)(2„i'2 +a'2) 



(33) 



Formally, ( po[ ) is an integral equation of the second kind. 
However, K, is not compact - the infinite integration range and 
the singularity on the a = a' diagonal prevent this -, and thus 
the solution theory for it is significantly different from what one 
has in the Fredholm case. We will return to this issue later. For 
now, let us just follow Zang in noting that uj factors out of the 
kernel, and thus the value of uj is irrelevant for the solvability 
of this equation. 

The special case of = has not been treated so far. It 
is interesting for two reasons: For one, when searching for the 
limit of axisymmetric stability, it suffices to consider uj — 0, 
since axisymmetric patterns cannot rotate anyway and consid- 
ering growth is not necessary for establishing the stability mar- 
gin, and for two, "neutral" modes (for which u; = even when 
TO 7^ 0) were the only modes the existence of which ERII and 
Zang report for the self-consistent disk. 

The expression ( p2[ ) cannot be taken to the limit a; ^ 0, 
because general has an essential singularity at = 

0, which according to Picard's theorem implies that the image 
of an arbitrarily small neighbourhood of is the entire complex 
plane with the possible exception of one point. Instead, one 
has to substitute uj — before the integration over L^. This 
simplifies the terms quite significantly, resulting in 



/Co(a') = 2J^(a') 



(34) 



The integral equation has now collapsed into the algebraic 
equation 



= /Co(a') - 1, 



(35) 



hia') 



which clearly simplifies solving the equation by at least an or- 
der of magnitude. 
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4. The cut-out disk 

The Mestel disk has three rather bizarre properties: Its mass is 
infinite, the surface density is not bounded in the centre, and it 
does not possess a natural time-scale. All these shortcomings 
can be cured by a relatively simple measure that was akeady 
employed by Zang. The idea is to immobilise a radially varying 
part of the mass, declaring it as, say, halo or bulge. Thus one 
keeps the simple structure of the disk, but cuts the singularity 
in the centre and possibly the long mass tail out to infinity as 
far as the Boltzmann equation is concerned. Also, the cut-out 
will break the self-similarity and thus introduces a time scale. 

As an added benefit, the ratio of active to total mass and 
thus (at radially constant v elocity dispersion) Toomre's sta- 
bility parameter Q (Toomre 1964 1 becomes a function of the 
radius, opening the possibility to introduce a Q-barrier in the 
disk. In such a Q-barrier, incoming waves are refracted away 
from the centre and thereby evade their absorption at the inner 
Lindblad resonance. In common se mi-ph enomenological theo- 
ries of spiral structure (Bertin et. al. 1989 ), this Q-barrier serves 
as a part of a feedback loop, to be clos ed by, e.g., the corotation 
at which swing amplification (Toomre 1981) or related mecha- 
nisms would not only provide reflection but even amplification 
of wave packets. 



In the shea ring sh eet (Goldreich & Lynden-Bell 1965 , Ju- 
lian & Toomre 1966 1, it is possible to show rigorously that the 
introduction of something quite similar to a Q-barrier leads to 
the emergence of stationary modes where no su ch mo des are 
present in the classic border-less analyses (Fuchs 2001 ). 

Technically, the immobilisation of the mass is done by mul- 
tiplying the distribution function with a cut-out function H{r), 
with < H{r) < 1, while keeping the Hamiltonian and thus 
the rotation curve unchanged. Since the main objective of our 
investigation of the cut-out disk was an assessment of the qual- 
ity of our approximation, we used the cut-out function used by 
ERI (who in turn basically used what Zang had suggested). In 
our case, the introduction of the cut-out has the effect of elim- 
inating the border terms in the integration, but only if one 
uses Zang's doubly cut-out disk (i.e., the cut-out function tends 
to zero for both r — > and r oo). Zang's family of cut- 
out functions is parametrised by two positive integers AI and 
TV, the cut-out indices, describing the steepness of the cut-outs 
in the outskirts and centre of the disk, respectively. Since both 
of these integers appear as the order of poles of the integrand 
in the integration, and the residues at these poles get very 
messy indeed as their order increases, we settled for the case 
M ^ N — 2, resulting in a cut-out function 



{LI 



Ll){Ll 



(36) 



where Lq and Lc are two parameters qualifying the locations 
of the inner and outer cut-outs; following ERII, we usually use 
Lq = 1 and Lc — 10. 

It should be noted that for further analytic exploration of 
the cut-out disk, this cut-out is not favourable, since it still 
leads to a very cumbersome kernel. A cut-out of H{Lz) — 



2X^X0/(^0 + ^z) results in a much more compact kernel and 
reproduces the overall behaviour the cut-out disk (due to its 
lower active surface density, the resulting disk has an even 
higher stability than one with Zang's cut-out, though). 

Since the cut-out function (|6|) only depends on Lz, most 
of the calculations done for the self-consistent disk just carry 
over The integration over becomes somewhat easier, since 
there are no border terms, but this advantage is more than offset 
by the sheer length of the terms and the occurrence of higher- 
order poles. The resulting integral equation has an "almost" 
compact kernel, with no singularity on the diagonal and a rel- 
atively sharp decline of the kernel as a and a' approach infin- 
ity. Consequently, ERI's naive approach of simply Fredholm- 
discretizing the integral equation and cutting off at some appro- 
priate a still works well. The only difference to ERTs method 
of solving this integral equation was that we found a spline- 
based quadrature to converge faster than the integration scheme 
proposed by them. 

Skipping the details - which can be found in Demleit- 
ner (2000 1 -, let us briefly compare the stability of disks de- 



scribed with our approximation with Read's findings. Since 
ERII mostly report results for singly cut-out disks which are 
not favourable for an anaylsis using our tools, we give the val- 



ues from Read's (1997) more comprehensive tables. Our disks 



become unstable to axisymmetric (m = 0) perturbations at 
Cmin = 0.274; Read only gives values for Lc = 100, where we 
find (Tniin = 0.349 versus Read's (Tmin — 0.327. This agree- 
ment to within roughly 5% is remarkable considering the rela- 
tively high temperature of this disk. The higher stability of the 
disk with lower Lc is mostly due to its lower active surface 
density. 

The one-armed (m — 1) modes investigated in ERII behave 
somewhat erratic in that at vanishing growth rate and 5Rtj 
in that there seem to be modes for almost arbitrary a. That, 
however, is a consequence of the fact that ERII performed 
these computations without an outer cut-out. Introducing an 
outer cut-out, the disk seems to become stable at a = 0.57, al- 
though clearly at that point the assumption of a razor-thin disk 
becomes quite unrealistic and Uc S> cr is severely violated. 

In the - at least historically - most interesting case of two- 
armed (m = 2) modes, our disks become unstable at aaun = 
0.214 with a pattern speed Q,p = uj /m = 0.548, whereas Read 
finds (Tmin = 0.205 and fip = a;/m = 0.542. Thus, not only 
is the qualitative result that the M = N = 2 cut-out disk is 
stable with respect to two-armed perturbations reproduced, the 
quantitative agreement again is far better than one might expect 
given the rather severe assumptions that go into our model. 



Read (1997 ) does not give stability limits for the m = 3 and 
m — A perturbations in doubly cut-out disks, but a comparison 
of the eigenvalue plots leads to the expectation that our results 
- fmin = 0.141 at J7p = 0.61 for m = 3 and cr,„i„ = 0.124 
at ftp = 0.588 for m = 4 - would again reproduce the results 
of the more exact calculations quite well. The m — 4 stability 
limit has to be taken with a grain of salt, since our formalism 
reproduces modes at high velocity dispersion as w ^ in ac- 
cordance with Read's findings. It is remarkable, though, that 
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these modes are not quenched by an outer cut-out Uke their 
m = brethren. 

The properties of the modes computed from the kernels 
are, by and large, identical to those discussed by ERII, and our 
model even shows more subtle features like the rapid overtak- 
ing of eigenvalues when m — 4. 

As to comparison to simulations, the fact that the expres- 
sions get very unwieldy for larger M and N is somewhat unfor- 
tunate. Sellwood & Evans (2001) point out that for the global 
mode to dominate over particle noise in simulations, it is highly 
desirable to maintain a high amplification, implying steep cut- 
outs and thus relatively large cut-out indices. The values cho- 
sen by Sellwood & Evans for their "standard" disk, = 4 
and M — 6, are beyond the reasonable range of applicabil- 
ity for our approach. However, they compute a doubly cut-out 
model that does not possess any global modes at Q = 1 (re- 
member that Q refers to the axisymmetric stability of the self- 
consistent disk) and state that the iV = 1 and M = 2 disk 
with Lq — 1 becomes stable with respect to to = 1 perturba- 
tions below Lc ~ 5. The rotation curve they use in this model 
is falling but only very slightly so, so that a disk with a com- 
pletely flat rotation curve should behave analogously. Indeed, 
ignoring the "rogue" eigenvalues (Demleitner 2000) - that be- 
come relevant in this model -, we find that even the to = 1 
mode gets quenched below Lc = 4.7. 

Summarising, our approximation does well in reproducing 
the results from the previous calculations for the cut-out disks. 
For now, the value of this insight is rather limited, since one 
still has to take recourse to numerical integration to solve the 
integral equation, and the numerics are computationally about 
as expensive as those in Zang's formalism, since the expression 
for the kernel is very large. Still, using a cut-out more suitable 
to our formalism results in more compact and manageable ex- 
pressions that might facilitate a closer look at the interplay be- 
tween the Q-barrier, the various resonances, and density waves 
in a full-disk analytical model. 

5. The self-consistent disk 

The axisymmetric stability of the self-consistent disk can be 
determined by solving the algebraic equation (DSb with a ker- 



be called rather peculiar: One-armed perturbations cannot be 
stabilised at all, two-armed modes set in below a = 0.155 (or 



cr = 0.149 according to Read (1997)), three-armed ones are 



nel still somewhat simpler than (34) that is the result of setting 
m = in (H) and yielding CTmin = 0.3860 with a = 
3.51, within 2% of Read's ( |1997| ) values of (T^in = 0.3780 and 
a = 3.46. 

Since the integral equation describing the self-consistent 
disk for non-axisymmetric modes is highly singular, Zang's 
formalism could not provide solutions for the full problem. 
However, it did yield a family of rather exotic modes charac- 
terised by a; = and consisting of a single logarithmic spiral 
- the neutral modes. It should be noted that these single log- 
arithmic spirals are scale-free, a property that sets them apart 
from from any other disturbance of the form (pi]). In our for- 
malism, they can be computed by solving the algebraic equa- 
tion (|5|) with the kernel (^4|). The resulting behaviour again 
agrees very well with the findings of ERII and Zang and has to 



extremely stable and do not set in down to at least a — 0.05, 
whereas four-armed neutral modes again cannot be stabilised 
at all. 

The neutral modes come in pairs, i.e., if a mode with log- 
arithmic wave number a exists, so does one with —a. This is 
a consequence of the anti-spiral theorem (Lynden-Bell & Os- 
triker 1967), since these neutral modes are not limits of grow- 
ing modes. The marginal modes in the cut-out disk and possibly 
existing modes with nonzero uj in the self-consistent disk, on 
the other hand, are limits of growing modes by the application 
of Landau's rule and thus not subject to the anti-spiral theorem. 

ERII, from their perspective of analysing the entire fam- 
ily of physically viable power-law disks, briefly hypothesised 
that the stability properties might be related to the emergence 
of closed orbits of a given symmetry; these would, according 
to her reasoning, absorb power from perturbations and stabilise 
the disks when present. In our model, no closed orbits at all ex- 
ist, since the ratio of epicyclic and orbital frequencies is fixed 
and irrational. As the model reproduces the stability properties 
found by ERII nevertheless, one has to look for a different ex- 
planation. 

Turning to (Q), one sees that the contribution from the 
terms except the one with the sum is fairly independent of 
TO, whereas the sum does significantly change its value with 
TO. The term to watch here is the TO'^/(m'^ — 2k^). When 
to' « V2k, the corresponding summand is much larger than 
its neighbours since the ratio Ik{x) and Ik+i{x) basically is k 
(plus some constant) and everything else does not depend on k; 
when m is not too large, that "resonant" summand will dom- 
inate the value of the entire series. The sign of this resonant 
summand is determined by the side relative to the nearest k on 
which to' — v^fc reaches its minimum - for the m = 3 se- 
ries, the near-resonance is at fc = 2 and thus one has a positive 
contribution with to' — ^/2k ~ 0.17 that enters into the kernel 
negatively (enhancing stability), whereas the resonance for the 
TO = 4 occurs for to' — \/2k = —0.24 and thus destabilises 
the disk. The next azimuthal harmonic that has an to close to a 
\/2k is 7, and it is again negative; not surprisingly, the to = 7 
harmonic is quite unstable as well. 

Unfortunately, k has no immediate meaning in physical 
terms, although (|l7]), suggests that the k-th terms of the se- 
ries in (Q) should determine the response in the k-th harmonic 
in the epicyclic angle wi . In the denominator of the resonance 
term to'^ — 2fc^, the factor \/2 is simply the epicyclic frequency, 
so that we seem to see a resonance between the frequency with 
which a body in the disk encounters a forcing, mil, and a har- 
monic of the epicyclic frequency, kn, which is a natural fre- 
quency a body wants to respond with. If the forcing is near res- 
onant with the response times the epicyclic frequency but has 
a slightly higher angular symmetry, the response is damped, 
whereas in the reverse case it is amplified, a behaviour common 
in forced oscillations. This resembles the situation at a Lind- 
blad resonance, with the difference that for one, the resonance 
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is not exact (if it were, the series would evaluate to infinity), 
and, for second, it occurs on the entire disk. 

Let us now turn to the general case of rotating modes. As 
has been hinted at above, the main discomfort with the notion 
of modes in the singular disk is that in its governing equation uj 
factors out, so that without further devices a disk is either stable 
or becomes unstable at all growth rates and pattern speeds at 
once. This seems weird enough to dismiss the possibility of 
such modes. On the other hand, this would imply that even a 
completely cold disk would be stable to all non-axisymmetric 
perturbations, which again would be quite an exotic situation. 



Quite recently, Goodman & Evans ( |1999| ) have tackled the 
problem employing the Jeans equations. The basic results of 
their work are that the Mestel disk's weirdness starts in its cen- 
tre, and once one fixes boundary conditions there, the Mes- 
tel disk becomes quite tame. Without that boundary condition, 
however, the problem is not well-posed. By the requirement 
that the centre does not absorb or emit energy, they curb down 
the two-dimensional continuum of solutions to a spectrum of 
one-dimensional continua, in which the ratio between growth 
rate and pattern speed is fixed and only the phase of uj remains 
unknown. Even this ambiguity can be removed by fixing the 
phase shift during the reflection of a wave at the disk's centre, 
thus reducing the spectrum to a countable set of discrete points; 
however, there is nothing in the equilibrium model that would 
allow the definition of this phase change. 

With a closed expression for the kernel of the Mestel disk, 
one is in a position to reconsider the stellar-dynamical equiva- 
lent of Goodman & Evans' problem using the tools of the the- 
ory of singular integral equations: Does the integral equation 
(30) admit continuous solutions, and if, for what disk parame- 
ters? 

The kernel ( ^2[ ) governing the self-consistent disk has a 
Cauchy-type singularity on the diagonal, implying that the op- 
erator / IC{a,a')^{a) da is not compact, so that the famil- 
iar Fredholm theory cannot be applied. For th e trea tment of 
integral equations of this type, Muskhelishvili ( 1953 ) and his 
coworkers have developed a theory. A nice summary is foun d 
in the more modern work of Polyanin & Manzhirov ( 1998 ). 
The following discussion is based on these works. 

The integral equation (BQ) has the form 







a(a')*(a') + — / $(a)da, (37) 

ITT Ji^ a — a 



where K is at least Holder continuous, L is the real axis (in- 
tegration over the real axis is henceforth implied for integral 
signs), and, of course, the integral is to be understood as a 
Cauchy principal value. In the present case, a{a') — T{a') — 1 
and K(q;, a') — «7r(a — Q;')/C(a, a'). The multiplication with 
[a — a') ensures that the kernel is bounded on the diago- 
nal, since the single zero in the denominator stemming from 
exp(27r(a — a')) — 1 now cancels out. The only part of the 
kernel that could jeopardise the condition of Holder continuity 
is the square root, but since only is present in the arguments 
of the square roots, K does indeed satisfy a Holder condition. 



It is convenient to rewrite ( |37| ) to the equivalent form 

= a{a')^a') + —b{a') f -^^da 
in I a ~ a' 



/ Kr(a, Q;')<I>(Q!)da, 

iTT J 



(38) 



where 



b{a) = K(a',a') and 

K(a,a')-K(a',a') 
Kr(Q:, a ) — — 



a — a' 

After this manipulation, is a compact operator, whereas 



(39) 



K$ := a(a')*(a') + —bia') 

ITT 



da — 0, 



(40) 



known as the dominant or characteristic part of the integral 
equation, decides whether or not the integral equation has a 
nontrivial solution. For our equation, b{a') — J- {a'). 

To find a solution of (p^, let us introduce an auxiUary func- 
tion 



^ ' 2m J 



{a-z) 



da, 



(41) 



defined for complex z. If $^(0;) is continuous - this is hence- 
forth assumed -, this function is piecewise analytic in the up- 
per and lower half planes but may be discontinuous across the 
real axis. The merit of this auxiliary function is that it allows 
a transformation of the problem of the solution of the integral 
equation ( p^ into a Riemann boundary value problem. 

This can be done by applying the Sokhotski-Plemelj formu- 
lae, that, in their formulation for the real axis, state that when 
$ and * are linked by ( |4l| ). 



^+{a') + m-{a') = — [ ^^da 
in J a — a' 



(42) 



hold for all a' e R. Here, ■^+{a') and are the limiting 

values of ^'(z) as z ^ a' from above and below the real axis, 
respectively. 

Inserting ( p^ into ( p^ and collecting the terms yields 



a{a') + b{a') 



(43) 



Going back to the integral equation for the singular disk, one 
has 



D{a')^-l/{2T{a')~l). 



(44) 



This represents a Riemann boundary value problem, asking (in 
the special case relevant here) for two functions and 'i>~ 
that are analytic on the upper and lower half planes, respec- 
tively, and that satisfy '^'+{a') = D{a')'^^ (a') for aU real a'. 

For reasons that will become clear below, we for the time 
being demand that D has neither zeroes nor poles on the real 
axis. Comparing (|5|) and (^4|), one sees that this will be the 
case exactly when neutral modes are admitted by the disk in 
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question, which strongly suggests that ( pO| ) has no solution in 
the regular case of a bounded and nonzero D, and at least one 
solution otherwise. With this concept mind, let us return to the 
analysis of the boundary value problem, for now keeping our 
attention on the case of a regular D. 

The D{a') treated here is the restriction of an analytic func- 
tion (in particular, by assumption it has no poles on the real 
axis), and thus is guaranteed to be Holder continuous. Further- 
more, -D(a') is a positive real function for real a'. This ensures 
that its index, the increment of its argument over the integration 
contour. 



Ind(i:>) = 



1 

2^ 



d'AigD{a), 



(45) 



is zero. This in turn implies that hi D is a well-defined (single- 
valued) function. Then one can take the logarithm of both sides 
of(g]), to arrive at 

In - Inl'-(a') = lni:>(a'). (46) 



Using the second of the Sokhotski-Plemelj formulae (| 
it follows that with 



G{z) 



1 

27ri 



\iiD{a) 



da 



(47) 



one can write down a solution to the boundary value problem 



(48) 



provided that the problem does have a solution. 

Thus one has D = /X^, and the boundary value prob- 
lem ( p3| ) can be written as 



X+{a') ^ X-{a')' 



(49) 



Now, the functions G by construction have no pole in their re- 
spective domains, and therefore X^ have no zeroes. If <i>^ are 
defined via an integral as in (^l]), this implies that they are ana- 
lytic in their respective half-planes. Thus, the functions on the 
two sides of (^9|) are analytic in their respective half-planes, 
and they are identical on the real axis. By Morera's theorem 
this implies that the right-hand side and the left hand side rep- 
resent the same analytic function, S{z), say. Now, what are the 
properties of S{z)l 

We already pointed out that S{z) is analytic on the en- 
tire complex plane. At infinity, one has ^^(oo) — ±^$(oo). 
Since we require any solution of the original integral equation 
(30) to vanish at infinity, one has 5'^(oo) = as well. Since 
D ^ 1 for large a, G approaches and X^ tends to 1. This 
implies that S vanishes at infinity. Since S is an entire func- 
tion, it follows from Liouville's theorem that S must be identi- 
cally zero. By the first Sokhotski-Plemelj formula one now has 
4>(a') = S'(a)(X++X~), and hence the homogenous integral 
equation has no nontrivial solution. This is a special case of the 
general result that is proven somewhat more rigorously in the 
literature cited above: A homogenous Cauchy integral equation 
with an index of zero is unsolvable. 



Loosening the conditions on D, we now allow it to have 
"poles" at some points Ui, which, for the sake of simplicity, 
are supposed to be poles of order one. The term "pole" has 
been put in quotes since it is not necessary that D is analytic 
in C \ {ai}, only that D{ai) behaves like (a' — ai)^^ at the 
exceptional points. 

Writing 



D'ia') 



n»=i(a - "i) 
the boundary condition ( p3| ) takes the form 

D'{a') 



*+(a') 



n:=i(a' 



(50) 



(51) 



where D' now satisfies the conditions we required of D in (p^]). 
Solving the boundary value problem for this D' and substitut- 
ing the result back into (|5l|), the equivalent of (^9|) now is 



X+{a') ~ X-(a')nr=i(«'-"0^ 



(52) 



where X and G are again defined by (^8|) and (^7|), except that 
now D is replaced by D' . Let us again call the analytic function 
defined in the left and right sides of ( p2| ) S{z). 

By the right side of (|5^, S now may have poles of order 1 
at ai but still has to vanish at infinity. Applying the generalised 
Liouville theorem, it follows that any solution of the integral 
equation must be of the form 



(53) 



Evidently, solutions of this type cannot vanish at oo, since such 
solutions will at best converge to X^ (oo) times a nonzero con- 
stant as a' ^ cx). 

Thus, provided that all singularities of S{z) are poles|], the 
self-consistent disk admits no growing modes at all. The ob- 
vious explanation for this bizarre absence of growing modes 
even in very cold disks is that growing modes would introduce 
scales into the scale-free disk. 

Stepping back to the mathematics for a moment, the fact 
that T is purely real on the real axis is crucial for that result, 
since even a minute complex contribution could provide for 
modes. For example, adding an ie to would not change the 
result that the index ( p3| ) of the problem ( |43| ) vanishes when 
T{a') < 1/2 for all a' , so that there still would be no modes 
above the stability limit of the neutral modes. Below that stabil- 
ity limit, however, the index would be nonzero and the integral 
equation would admit nontrivial solutions. 

Such an imaginary contribution would necessitate that the 
real part of the kernel ( ^2[ ) shows a singularity on the diagonal. 
We cannot offer an explanation how this might be linked to the 
scale-freeness of the disk at this point. 



Thi s assu mption is hard to prove, but generally true, which is why 
Pipkin ( 1991 ) calls it Gentlemen's Theorem Number 1 - a theorem no 



gentleman would question. 
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6. Conclusions 

In this paper we have derived an analytic expression for the 
kernel of the Mestel disk. In the cut-out disk, it reproduces the 
behaviour described by Zang and ERll both qualitatively and 
quantitatively surprisingly well. It turns out that the expressions 
for the kernel are quite large in the cut-out case, at least when 
one uses the classic cut-out functions suggested by Zang, so 
the computational effort involved in finding modes is about as 
large with our formalism as it is with the older seminumerical 
scheme. On the other hand, it seems that with carefully selected 
cut-out functions, the kernel could be simple enough for further 
analytic work. 

The kernel governing modes in the self-consistent disk is 
relatively compact. It allowed us to interpret the stability be- 
haviour of neutral modes in terms of global resonances be- 
tween the orbital (which coincides with the frequency of excita- 
tion in this case) and epicyclic frequencies. In the longstanding 
question of rotating modes we could establish that no rotating 
modes can exist above the stability limit of the neutral modes. It 
might seem surprising that rotating modes should know about 
the stability limit of the rotating modes, considering that the ra- 
tio of orbital and epicyclic frequency is no longer independent 
of the radius for rotating modes. However, due to the lack of a 
length scale in the Mestel disk, any mode is arbitrarily close to 
a non-rotating one. 

Even below the stability limit of the neutral modes, no ro- 
tating or growing modes exist in the perfectly scale-free disk 
under reasonable assumptions like continuity of the coefficient 
function in the expansion of the potential perturbation (pi]). 
This bizarre property is probably best explained by an inabil- 
ity of the disk to break its self-similarity and thus by a com- 
pletely artificial feature. In conclusion, it seems to us that the 
self-consistent disk's peculiar properties make its investigation 
something of a pedagogical (Goodman & Evans 1999) exer- 
cise, whereas the cut-out disk might yet prove to be a valuable 
testbed for wave dynamics in a fairly realistic full-disk setting. 
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